#############################################################################################
########### calculate exposure index for each species ########################################
#############################################################################################

#load packages:
library(R.utils)
library(raster)
library(sp)
library(sf)
library(stringr)
library(circular)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
MUOutDir<-paste(MyDir, "/Outputs2_Crossval/permutation.importance", sep="")
IndSppOutDir<-paste(MyDir, "/Outputs3_Final/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define 'for' variables:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
timecats<-c("Current","2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
#timecats<-c("Current","2050_Median","2090_Median")
#timecats<-c("Current")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-DatSum[DatSum$Listing_comment!="not listed",]
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]
#spSTART<-which (spp=="Micrurus_mertensi")
DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species

#################################################################################
#################################################################################
# plot trends (individual species lines, and overall line)

###################################
# 0. Individual Species Plots #####
###################################


myX<-c(1985,2050,2090)
#####################################################################
 DatSumMod$RangesizeSQKM_Current[DatSumMod$RangesizeSQKM_Current==0]<-0.01
 DatSumMod$SuitabilitySum_Current[DatSumMod$SuitabilitySum_Current==0]<-0.01
 DatSumMod$Exposure_Current[DatSumMod$Exposure_Current==0]<-0.01
 DatSumMod$Impact_Current[DatSumMod$Impact_Current==0]<-0.01

for (col in c(13:26,41:61)){
  DatSumMod[,col]<-  ifelse(DatSumMod[,col]==0,1,DatSumMod[,col])
}
#####################################################################
#for (sp in c("Bungarus_multicinctus","Naja_nigricollis","Agkistrodon_piscivorus","Agkistrodon_contortrix")){
   
for (sp in spp){
  
  print(sp)
  
  n<-which(spp==sp)
  
  #1. make circular plot:
  
  png(filename=paste(PlotDir,"/",sp,"_CCShift.png",sep=""), units="in", width=10, height=10, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(1,1), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  DirDeg<-    c(DatSumMod$DirDeg_2050_Q10[n],
                DatSumMod$DirDeg_2050_Median[n],
                DatSumMod$DirDeg_2050_Q90[n],
                DatSumMod$DirDeg_2090_Q10[n],
                DatSumMod$DirDeg_2090_Median[n],
                DatSumMod$DirDeg_2090_Q90[n])
  DirDist<-   c(DatSumMod$DistM_2050_Q10[n],
                DatSumMod$DistM_2050_Median[n],
                DatSumMod$DistM_2050_Q90[n],
                DatSumMod$DistM_2090_Q10[n],
                DatSumMod$DistM_2090_Median[n],
                DatSumMod$DistM_2090_Q90[n])/1000
  DirDegCirc<- circular(DirDeg, type = "directions", 
                        units = "degrees",
                        modulo = "asis", 
                        rotation = "clock")
  
  
  #make base plot with package plotrix
  radial.plot(
    lengths=c(DirDist), 
    radial.pos=c(DirDeg*(pi/180)), 
    label.pos=c(0,pi/2,pi,pi+pi/2),
    labels=c("North","East","South","West"),
    lwd=c(2,4,2,2,4,2),
    lty=c(3,1,5,3,1,5),
    start=(pi/2),
    clockwise=TRUE,
    line.col=colors,
    radial.lim = c(0,max(DirDist)), #range of grid circle
    main=paste(gsub("_"," ",sp), "Predicted Range Shift"),
    show.centroid=FALSE,
    show.grid.label=1, #put the concentric circle labels going down
    show.radial.grid=TRUE,
    rp.type="n",
    grid.unit=" km"
  )
  
  radial.plot(add=TRUE,
              lengths=c(0,DirDist[1:3]), 
              radial.pos=c(0,DirDeg[1:3]*(pi/180)), 
              start=(pi/2),
              clockwise=TRUE,
              rp.type="p",
              line.col=rgb(col2rgb("white")[1],col2rgb("white")[2],col2rgb("white")[3],max=255,alpha=0),
              poly.col=rgb(col2rgb("grey60")[1],col2rgb("grey60")[2],col2rgb("grey60")[3],max=255,alpha=75),
              radial.lim = c(0,max(DirDist))) #range of grid circle
  
  radial.plot(add=TRUE,
              lengths=c(0,DirDist[4:6]), 
              radial.pos=c(0,DirDeg[4:6]*(pi/180)), 
              start=(pi/2),
              clockwise=TRUE,
              rp.type="p",
              line.col=rgb(col2rgb("white")[1],col2rgb("white")[2],col2rgb("white")[3],max=255,alpha=0),
              poly.col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75),
              radial.lim = c(0,max(DirDist))) #range of grid circle
  
  # plot.circular(DirDegCirc, pch = 16, cex = 1, stack = TRUE,
  #               axes = TRUE, start.sep=0.05, sep = 0.025, shrink = 1,
  #               bins = 360, ticks = TRUE, tcl = 0.025, tcl.text = 0.125,
  #               col = "red", tol = 0.04, uin = NULL,
  #               xlim = c(-max(DirDist), max(DirDist)), ylim = c(-max(DirDist), max(DirDist)), digits = 2, units = NULL,
  #               template = NULL, zero = (pi/2), rotation = NULL, 
  #               main = NULL, sub=NULL, xlab = "", ylab = "", 
  #               control.circle=circle.control(type="l",bg="lightgrey",col="grey20",lwd=2))
  # 
  
  #Draw arrows with package circular:
  arrows.circular(DirDegCirc,DirDist, x0 = 0, y0 = 0, na.rm = FALSE,
                  plot.info = NULL, zero = (pi/2),col=colors,lwd=c(2,4,2,2,4,2),
                  lty=c(3,1,2,3,1,2))
  
  radial.plot(add=TRUE,
              lengths=c(DirDist+(max(DirDist,na.rm=TRUE)/20)), 
              radial.pos=c(DirDeg*(pi/180)), 
              start=(pi/2),
              clockwise=TRUE,
              radial.lim = c(0,max(DirDist)), #range of grid circle
              rp.type="t",
              boxed.radial=FALSE,
              point.symbols=c("2050\nQ10","2050\nMedian","2050\nQ90","2090\nQ10","2090\nMedian","2090\nQ90"))
  
  dev.off()
  
  #2. make species line graphs
  
  RangeMed<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossMed<-c(1,
              1-DatSumMod$Loss_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainMed<-c(1,
              1+DatSumMod$Gain_2050_Median[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Median[n]/DatSumMod$RangesizeSQKM_Current[n])
  
  RangeQ10<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossQ10<-c(1,
              1-DatSumMod$Loss_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainQ10<-c(1,
              1+DatSumMod$Gain_2050_Q10[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Q10[n]/DatSumMod$RangesizeSQKM_Current[n])
  
  RangeQ90<-c(DatSumMod$RangesizeSQKM_Current[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              DatSumMod$RangesizeSQKM_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  RLossQ90<-c(1,
              1-DatSumMod$Loss_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              1-DatSumMod$Loss_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  RGainQ90<-c(1,
              1+DatSumMod$Gain_2050_Q90[n]/DatSumMod$RangesizeSQKM_Current[n],
              1+DatSumMod$Gain_2090_Q90[n]/DatSumMod$RangesizeSQKM_Current[n])
  #####################################################################
  SuitMed<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Median[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Median[n]/DatSumMod$SuitabilitySum_Current[n])
  SuitQ10<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Q10[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Q10[n]/DatSumMod$SuitabilitySum_Current[n])
  SuitQ90<-c(DatSumMod$SuitabilitySum_Current[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2050_Q90[n]/DatSumMod$SuitabilitySum_Current[n],
             DatSumMod$SuitabilitySum_2090_Q90[n]/DatSumMod$SuitabilitySum_Current[n])
  #######################################################################
  RiskMed<-c(DatSumMod$Risk_Current[n]/DatSumMod$Risk_Current[n],
            DatSumMod$Risk_2050_Median[n]/DatSumMod$Risk_Current[n],
            DatSumMod$Risk_2090_Median[n]/DatSumMod$Risk_Current[n])
  #######################################################################
  ExpMed<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Median[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Median[n]/DatSumMod$Exposure_Current[n])
  ExpQ10<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Q10[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Q10[n]/DatSumMod$Exposure_Current[n])
  ExpQ90<-c(DatSumMod$Exposure_Current[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2050_Q90[n]/DatSumMod$Exposure_Current[n],
            DatSumMod$Exposure_2090_Q90[n]/DatSumMod$Exposure_Current[n])
  #######################################################################
  ImpMed<-c(DatSumMod$Impact_Current[n]/DatSumMod$Impact_Current[n],
            DatSumMod$Impact_2050_Median[n]/DatSumMod$Impact_Current[n],
            DatSumMod$Impact_2090_Median[n]/DatSumMod$Impact_Current[n])

  
  
  ###########################################
  ###########################################
  ###########################################
 
  png(filename=paste(PlotDir,"/",sp,"_CCTrends.png",sep=""), units="in", width=14, height=7, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(1,2), mar = c(5,5,2,2), oma=c(0.2,0.2,0.2,0.2), mgp=c(3.2,1.4,0))
  # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  plot(myX,RangeMed, ylim=c(0,max(ExpQ90,SuitQ90,2,na.rm=TRUE)),xlim=c(1980,2090),col="white",bg="white", 
       bty="l", pch=21, lwd= 4, ylab="fraction of current", xlab="year",lty=1, main=NULL,xaxt="n",yaxt="n",cex.lab=2)
  
  axis(side=1, at=c(2010,2050,2090), lwd=1, cex.axis=2)
  axis(side=2, at=c(0,1,2,3,4), lwd=1, cex.axis=2)
  
  #plot zero line for reference
  lines (c(1980,2050,2090),c(1,1,1), col="grey", lwd= 2, lty=1)
  
  #2.a. plot % change in range size (black, solid thin)
  polygon(c(myX,rev(myX)),c(RangeQ10,rev(RangeQ90)), border = NA,
          col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))
  #points(myX,RangeMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,RangeMed, col="black", lwd= 4, lty=1)
  lines (myX,RangeQ10, col="black", lwd= 2, lty=2)
  lines (myX,RangeQ90, col="black", lwd= 2, lty=2)
  
  #2.b. plot % of range lost (1985=0) (blue, solid thin)
  polygon(c(myX,rev(myX)),c(RLossQ10,rev(RLossQ90)), border = NA,
          col=rgb(col2rgb("royalblue")[1],col2rgb("royalblue")[2],col2rgb("royalblue")[3],max=255,alpha=75))
  #points(myX,RLossMed, col = "black", bg="royalblue",  pch=21, lwd= 1)
  lines (myX,RLossMed, col="royalblue", lwd= 4, lty=1)
  lines (myX,RLossQ10, col="royalblue", lwd= 2, lty=2)
  lines (myX,RLossQ90, col="royalblue", lwd= 2, lty=2)
  
  #2.c. plot % of range gained (1985=0) (red, solid thin)
  polygon(c(myX,rev(myX)),c(RGainQ10,rev(RGainQ90)),border = NA, 
          col=rgb(col2rgb("orangered")[1],col2rgb("orangered")[2],col2rgb("orangered")[3],max=255,alpha=75))
  #points(myX,RGainMed, col = "black", bg="orangered",  pch=21, lwd= 1)
  lines (myX,RGainMed, col="orangered", lwd= 4, lty=1)
  lines (myX,RGainQ10, col="orangered", lwd= 2, lty=2)
  lines (myX,RGainQ90, col="orangered", lwd= 2, lty=2)
  
  
  plot(myX,RangeMed, ylim=c(0,max(ExpQ90,SuitQ90,2,na.rm=TRUE)),xlim=c(1980,2090),col="white",bg="white", 
       bty="l", pch=21, lwd= 4, ylab="", xlab="year",lty=1, main=NULL,xaxt="n",yaxt="n",cex.lab=2)
  
  axis(side=1, at=c(2010,2050,2090), lwd=1, cex.axis=2)
  axis(side=2, at=c(0,1,2,3,4), lwd=1, cex.axis=2)
  
  
  #plot zero line for reference
  lines (c(1970,2050,2100),c(1,1,1), col="grey", lwd= 2, lty=1)
  
  #2.d. plot % change in suitability/abundance (dashed)
  polygon(c(myX,rev(myX)),c(SuitQ10,rev(SuitQ90)), border = NA,
          col=rgb(col2rgb("darkgreen")[1],col2rgb("darkgreen")[2],col2rgb("darkgreen")[3],max=255,alpha=75))
  #points(myX,SuitMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,SuitMed, col="darkgreen", lwd= 4, lty=1)
  lines (myX,SuitQ10, col="darkgreen", lwd= 2, lty=2)
  lines (myX,SuitQ90, col="darkgreen", lwd= 2, lty=2)
  #2.e. plot % of abundance lost (1985=0)?
  #2.f. plot % of abundance gained (1985=0)?
  
  #2.g. plot % change in human exposure (dotted)
  polygon(c(myX,rev(myX)),c(ExpQ10,rev(ExpQ90)), border = NA,
          col=rgb(col2rgb("orange")[1],col2rgb("orange")[2],col2rgb("orange")[3],max=255,alpha=75))
  #points(myX,ExpMed, col = "black", bg="black",  pch=21, lwd= 1)
  lines (myX,ExpMed, col="orange", lwd= 4, lty=1)
  lines (myX,ExpQ10, col="orange", lwd= 2, lty=2)
  lines (myX,ExpQ90, col="orange", lwd= 2, lty=2)
  #2.h. plot % new exposure (1985=0)?
  #2.i. plot % lost exposure (1985=0)?
  
  # #2.j. plot % change in impact (solid thick)
  # points(myX,ImpMed, col = "black", bg="black",  pch=23, lwd= 1)
  # lines (myX,ImpMed, col="black", lwd= 2, lty=1)
  # #2.k. plot % new impact (1985=0)
  # #2.l. plot % lost impact (1985=0)
  
  #plot risk
  
  #mtext(paste(gsub("_"," ",sp),sep=""), font= 4, side = 3, line = 1, outer = TRUE, cex= 1.5)
  
  dev.off()
  
  }


###############################
# 3. AUC ######################
###############################

{
  png(filename=paste(PlotDir,"/AUCs.png",sep=""), units="in", width=10, height=10, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(1,1), mar = c(1,1,3,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
  
  plot(log(DatSumMod$MUrecs+1), DatSumMod$TrainingAUC, ylim=c(0,1),xlim=c(0,max(log(DatSumMod$MUrecs+1))),col="black",bg="white", 
       type="p", bty="n", axes=TRUE, pch=21, lwd= 2, main="Model AUCs",cex=1.5)
  
  points(log(DatSumMod$MUrecs+1),DatSumMod$TestAUC, col="black",     pch=8, lwd= 2)
  
  dev.off()
}
 #############################
##############################
##############################
